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Abstract 

Steady state plastic flows have been compared to developed turbulence be- 
cause the two phenomena share the inherent complexity of particle trajectories, 
the scale free spatial patterns and the power law statistics of fluctuations. The ori- 
gin of the apparently chaotic and at the same time highly correlated microscopic 
response in plasticity remains hidden behind conventional engineering models 
which are based on smooth fitting functions. To regain access to fluctuations, 
we study in this paper a minimal mesoscopic model whose goal is to elucidate 
the origin of scale free behavior in plasticity. We limit our description to fee type 
crystals and leave out both temperature and rate effects. We provide simple il- 
lustrations of the fact that complexity in rate independent athermal plastic flows 
is due to marginal stability of the underlying elastic system. Our conclusions are 
based on a reduction of an over-damped visco-elasticity problem for a system 
with a rugged elastic energy landscape to an integer valued automaton. We start 
with an overdamped one dimensional model and show that it reproduces the main 
macroscopic phenomenology of rate independent plastic behavior but falls short 
of generating self similar structure of fluctuations. We then provide evidence that 
a two dimensional model is already adequate for describing power law statistics of 
avalanches and fractal character of dislocation patterning. In addition to capturing 
experimentally measured critical exponents, the proposed minimal model shows 
finite size scaling collapse and generates realistic shape functions in the scaling 
laws. 

Keywords: plasticity, dislocations, self-organized criticality, Frenkel-Kontorova 
model, statistics of avalanches, intermittency, power laws 
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1. Introduction 

It is known that crystalline and amorphous solids flow plastically when macro- 
scopic stresses exceed certain thresholds. In plasticity, thermal relaxation may of- 
ten be neglected while the driving can still be viewed as quasi-static. In this case 
the flow is an athermal process with rate independent dissipation. The singularity 
of this dissipative mechanism distinguishes plasticity from the regular dissipative 
phenomena such as viscosity or heat conductivity, and places it instead in the class 
of dry friction and fracture [1, 2, 3, 4, 5]. 

At the macroscale rate independent plasticity can be perceived as a smooth 
process, whereas at microscale the plastic flow is known to exhibit fluctuations 
revealing isolated dissipative events [7, 8]. In crystals these fluctuations can be 
linked to nucleation and collective depinning of dislocations interacting through 
long range elastic stress fields [9, 10]. In amorphous and random granular systems 
the nature of the 'quantum' of plastic strain is still debated [11, 12, 13, 14, 15] and 
therefore in what follows we mostly limit our discussion to crystals. 

The smooth description of crystal plasticity adopted in engineering theories 
presumes spatial homogenization and time averaging. In some cases (bcc metals, 
tetrahedral covalent crystals, etc.) the obstacles are strong, the dislocation inter- 
action is weak and the plastic flow can be viewed as a sum of uncorrelated events. 
The collective effects can then be neglected making classical homogenization an 
appropriate tool. Here we are interested in an alternative case (fee metals, hep 
crystals with basal glide, etc.) when hardening can be neglected, dislocational mo- 
bility is high and the elastic interaction among distant dislocations is significant. 
Then a specific collective behavior at the macroscale emerges as a manifestation of 
many correlated events at the microscale [16, 17, 18, 19]. In particular, such plas- 
tic flows exhibit in the steady state irregular isolated bursts and reveal apparently 
random localized active slip volumes, with both spatial and temporal fluctuations 
spanning all scales. The temporal intermittency manifests itself through acoustic 
emission with power law statistics of avalanches [20, 21, 22, 23] while spatial 
fluctuations take the form of fractal dislocational cell structures [24, 25, 26]. 

The ensuing macroscopic dynamics is usually characterized as critical because 
of the reasons that will become more clear in what follows. In physical terms the 
resulting behavior can be described as glassy. Since dislocational motion proceeds 
in a self similar manner in both space and time the importance of such notions as 
average dislocation density and average velocity in such regimes become ques- 
tionable. Those notions, however, remain relevant in non-critical plastic regimes 
where dislocational microstructures are regular and are characterized by particular 
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scales which depend on the level of loading, e.g. [27]. 

Despite the over-schematic description of dislocational activity in classical 
engineering theory, it has been remarkably successful in capturing the most im- 
portant plasticity phenomenology such as yield stress, work hardening and shake- 
down [29, 30]. However, the fractal patterning and the peculiar scale free structure 
of the temporal fluctuations remain invisible at the macroscale and are missed by 
the approaches based on 'intuitive' spatial and temporal homogenization. As a 
result, a quantitative link between the fitting functions used in phenomenologi- 
cal plasticity and the microscopic picture of a crystalline lattice with moving line 
defects remains elusive. This gap has been a major obstacle on the way of quanti- 
tative estimate of plastic response for artificially designed materials. 

The emergence of power laws suggests that the relation between microscopic 
and macroscopic pictures of plastic flow is rather complex and more akin to tur- 
bulence [31, 18, 32] than to elasticity where classical homogenization of lattice 
models is usually sufficient to predict macroscopic response [33, 34]. It is also 
clear that powerful methods of equilibrium statistical mechanics implying ergod- 
icity and absence of correlations are not applicable for the description of small 
scale plasticity which appear to be highly nonequilibrium and strongly correlated 
phenomenon. The main problem with applying the classical continuum approach 
is that spatial renormalization symmetry (self- similarity) excludes separation of 
scales while non gaussian statistics prevents local description. More generally, 
we still lack analytic tools allowing quantitative description of evolutionary phe- 
nomena that depend on rare events spanning the whole domain and have con- 
vergence problems handling the uncertainty which grows algebraically with time 
which qualifies such systems as located at the 'edge of chaos' [28]. 

The fine structure of long range correlations in plasticity has been extensively 
studied in recent experiments which unambiguously established the power law 
statistics of avalanches and produced numerous examples of the scale-free nature 
of dislocation patterns [20, 21, 22, 23, 35, 36, 37, 38, 16, 39, 40, 41, 42]. Similar 
phenomena of self organization towards criticality have been also observed in 
other mechanical systems operating in nonequilibrium steady regimes including 
friction, fracture, porous and granular flows and martensitic transformations [43, 
44, 45, 46, 47]. 

While the general mechanisms for the formation of scale free correlations re- 
main a major subject of research [48], several important ingredients of dynamic 
criticality, such as marginal stability, quasi-static driving, threshold type nonlin- 
earity and nonlocal feedback, have been identified [49, 50]. All these compo- 
nents are present in crystal plasticity. For instance, marginal stability reveals itself 
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through the fact that infinitesimal perturbations not only trigger small scale rear- 
rangements of the ensemble of locked dislocations but can also induce a global 
transition between the jammed and the flowing regimes. The long range interac- 
tions responsible for the feedback have elastic nature, suggesting that, somewhat 
paradoxically, elasticity plays a fundamental role even in 'ideal' plasticity. 

The experimental findings of plastic criticality in crystals have been supported 
by a number of numerical simulations reviewed in [51, 52]. The two main ap- 
proaches are: the discrete dislocation dynamics (DDD) accounting for disloca- 
tion interaction on different slip planes [53, 54, 53, 55, 56] and the pinning- 
depinning models considering dislocations dynamics on a single slip plane [57, 
58, 59]. Discrete mesoscopic models of crystal and amorphous plasticity implying 
some averaging and providing an effective description of the pinning-depinning 
and jamming-unjamming processes have also been shown to generate power law 
statistics of avalanches with realistic exponents [62, 63, 64, 66, 67, 15]. Among 
those we particular emphasize the models of earthquakes aimed at reproducing 
the Gutenber-Richter law because the rapture-healing phenomena on a preexist- 
ing fault are very similar to 2D plasticity showing almost the same critical ex- 
ponents [65, 60, 61]. In view of the implied analogy with turbulence, of con- 
siderable interest are also the models dealing with dynamics of continuously dis- 
tributed dislocations and capturing the scale free effects already in a PDE setting 
[69, 70,71,72, 73,32, 74]. 

Despite their ability to reproduce experimental findings, the above models can- 
not provide theoretical insight into the origin of criticality because they lack an- 
alytical transparency and depend on large scale numerical experiments. On the 
other hand, the physical adequacy of the proposed numerical schemes may still 
be disputed in some basic details. For instance, in case of DDD models, the fast 
microscopic topological changes (nucleation, annihilation, interaction with ob- 
stacles, etc) are treated by auxiliary hypotheses formulated phenomenologically 
in terms of local stresses without addressing the mechanism of barrier crossing 
[75, 77]. Likewise, the theory of continuously distributed dislocations have diffi- 
culty representing kinks, locks, jogs and other individual entanglements with non- 
trivial topology. The numerical models can be made more adequate by involving 
various quasi-continuum schemes, however, it is clear that in order to understand 
plastic criticality mathematically, the existing models have to be simplified rather 
than further complicated. 

The anticipation of a simple and transparent theory is based on the idea that 
truly scale free behavior is independent of either microscopic or macroscopic de- 
tails of the system. The statistics of dislocation avalanches should then be an 
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intrinsic feature of a particular crystal class described in terms of symmetry and 
dimensionality and not affected by specific material parameters and details of the 
loading geometry. It is then of interest to formulate simple prototypical models 
representing particular universality classes of plastic flows. The simplicity is un- 
derstood here in the sense that the model is amenable to rigorous mathematical 
study while still capturing not only the observed exponents but also the character- 
istic shape functions in scaling relations [49, 78]. 

At present, the only analytically tractable approaches to criticality in driven 
distributed systems are the models based on branching processes [44], the mean 
field model [79], the elastic depinning model amenable to FRG methods [80] 
and the Abelian sand pile type automata with integer valued fields [81, 82]. The 
challenge is to reduce a realistic model of plasticity to one of these analytically 
transparent types. 

Since the relation between the prototypical branching processes and the dy- 
namical phenomena in distributed systems remains rather imprecise, we will not 
pursue this line of modeling here. In contrast, the observed critical exponents 
for crystal plasticity are so close to the mean field values that claims have been 
made that 3D plasticity is in the mean field universality class. This hypothesis is 
supported by the estimates of the critical dimension in various systems with long 
range elastic interactions [49, 83, 84, 85]. However, the question of validity of 
this approximation remains open because the match of the measured exponents is 
not perfect and the observed avalanche shapes appear to be less symmetric than 
the theoretical predictions [86]. Also, the mean field values of exponents have not 
been supported by numerical experiments with colloidal crystals [87] and within 
phase field crystal model [68]. 

The depinning models have a considerable relevance for plasticity because 
they study quasi- statically driven elastic objects interacting with a set of randomly 
distributed obstacles. The corresponding theoretical critical exponents for the case 
of dislocations have been computed with high precision [88, 89, 51, 90] and con- 
firmed by molecular dynamics studies [91]. These models, however, may not be 
in the same universality class as crystal plasticity because the critical exponents 
are different from those observed in plasticity experiments. Most probably, the 
origin of divergence is the neglect in depinning models of important interactions 
among dislocations on different slip planes. 

Left with the only remaining option, we study in this paper the possibility to 
reduce the conventional dislocation mediated plasticity model to an integer val- 
ued Abelian automaton. The Abelian symmetry means that the order, in which 
the eligible 'spin variables' are updated during the avalanche, is irrelevant for the 
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avalanche outcome. Some Abelian automata are amenable to rigorous analysis 
due to partial ergodicity on a set of recurrent sequences which were fully charac- 
terized in several cases [81, 82, 92, 93, 94]. The aim of this paper is to show that 
a reduction of a multi-slip-plane dislocational dynamics to an integer automaton 
with some form of Abelian symmetry allows one to reproduce the experimentally 
observed statistics of avalanches. An analytical study of the ensuing automaton 
will be reported elsewhere. 

2. Summary of the main results 

In search for the simplest representation of plasticity universality class(es) we 
begin with a series of one dimensional models allowing one to reproduce rate 
independent dissipation and hysteresis. Then we move to two dimensional models 
capturing also intermittency and power law statistics. 

Our starting point is an assumption that the micro-scale dynamics is over- 
damped and that the rate of loading is much slower than the rate of viscous relax- 
ation but much faster than the rate of thermal relaxation. Such models belong to 
the class of AQS (athermal quasi-static) systems that have been found relevant for 
the description of many rate independent phenomena from wetting to magnetism 
[95, 96, 97, 98, 99, 100]. 

The simplest ID lattice model of this type, representing a chain of interact- 
ing particles (or a deck of interacting rigid cards), can be developed for trans- 
formational plasticity of shape memory alloys. To obtain in this setting a plas- 
tic response at the macro-level we must assume that the interaction potential has 
a double well structure which fully characterizes the local 'mechanism' of the 
transformation. The presence of local bi-stability, which is also relevant for the 
description of the 'easy glide' stage in crystal plasticity, places such mechanical 
system into the class of snap-spring lattices. The non-convexity of local interac- 
tions may be equally relevant for amorphous plasticity operating with a concept of 
STZ (shear transformation zones), which implies a multi-well structure of an ef- 
fective potential [5, 12]. We emphasize that the discreteness at this stage has to be 
understand broadly as describing a generic structural inhomogeneity, for instance, 
a grain structure. 

In the simplest chain with nearest neighbor (NN) interactions the elastic el- 
ements are subjected to common stress field. The resulting mean field model 
captures the basic mechanism of rate independent dissipation and can be used to 
illustrate the idea of marginal stability of the underlying elastic system which is 
ultimately responsible for the plastic yield [8, 101]. A more realistic model must 
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incorporate short range interactions and the simplest extension in this direction is 
a chain with next to nearest neighbors (NNN). The latter may be of ferromagnetic 
or anti-ferromagnetic type [102, 103] and in both cases the resulting mechanical 
system can be viewed as a soft spin generalization of the classical Ising model. 

In the case of ferromagnetic interactions penalizing gradients we show that 
that the NNN model captures the important difference between the nucleation and 
the propagation thresholds and describes the characteristic nucleation peak [104]. 
In the case of anti-ferromagnetic interactions, favoring lattice scale oscillations, 
such model can reproduce periodically distributed shear bands. We find, however, 
that even in the presence of gaussian quenched disorder, the internal dynamics 
in all these ID models is still too simple and none of them can generate scale 
free correlations in the overdamped setting. The situation does not improve if the 
double-well NN potential is replaced by a periodic potential. 

Despite their failure to capture criticality, the ID models are very important 
because they allow one to show that in the limit of quasi-static driving a viscous 
evolution in a rugged energy landscape takes the form of a stick slip dynamics 
adequately described by a discrete automaton. The replacement of the fast stages 
of dynamics by jumps leads to considerable simplification of the computational 
problem. It also shows that in quasi-static setting all dissipative mechanisms pro- 
ducing the same set of jumps are equivalent. 

To simplify the system even further we replace a smooth multi-well potential 
by its piece-wise parabolic analog allowing one to solve the elastic problem ana- 
lytically when the 'phase configuration' is known. Since the phase configuration 
is prescribed by an integer valued field, such 'condensation' of elastic problem 
allows one to reduce the original smooth dynamical system to an integer val- 
ued automaton of a sand pile (threshold) type. Because of the long range nature 
of elastic interactions, the discharge rules in this automaton are nonlocal which 
makes its rigorous analysis challenging. 

While minimizing out the elastic fields in plasticity problems is a rather com- 
mon approach leading to a reduced description in terms of either plastic distortion 
or dislocation density [105, 106, 76, 71, 32, 107], the proposed integer automaton 
representation appear to be new. Its salient feature is that in place of straightfor- 
ward time discretization of continuum dynamics [108] we utilize inherent tempo- 
ral discreteness hidden behind the conventional gradient flow dynamics (see also 
[109, 8, 101]). 

The next level of complexity is exemplified by 2D lattice models. Here we 
abandon the double-well framework with its prototypical defects and consider 
the simplest setting allowing one to represent actual dislocations which interact 
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through realistic long range elastic fields and are capable of generating collective 
effects. As it has been repeatedly proposed in the literature, the two dimensional 
models of plasticity may be sufficient to capture the corresponding universality 
class at least for FCC and hexagonal crystals [53, 63, 112]. For similar models of 
criticality associated with transformational plasticity, see [109, 110, 111]). 

In general terms, our 2D model represents a cross section of a crystal with a 
single slip geometry. More precisely, we consider a set of parallel edge disloca- 
tions which can move only in horizontal direction. This implies that the crystal 
is highly anisotropic with the deformation in one direction strongly constrained. 
As a partial justification we mention that intermittency has been mostly observed 
under single slip conditions [63, 112]. 

The ensuing discrete model with continuous dynamics can be viewed as a par- 
allel set of coupled overdamped FK chains [113, 1 14, 115, 116, 1 17]. It has been 
shown before that this model allows for generation and annihilation of dislocations 
and that it describes adequately their long range elastic interactions [118, 119]. 
Most importantly, both the propagation and the nucleation/anihillation mecha- 
nisms in this model are associated with reaching the same strain thresholds. In 
particular, this model correctly accounts for Peierls stress associated with under- 
lying periodicity, describes adequately sufficiently slow dislocational kinetics and 
does not require special rules governing the interaction of the dislocations with 
quenched disorder [117, 118]. At the same time, this model still neglects tensorial 
nature of elasticity and ignores such important physical effects as cross-slip and 
climbing which may affect critical exponents. 

In our numerical simulations we assume periodic boundary conditions and 
submit the 2D lattice to quasi-static cycling in the hard device by imposing pe- 
riodic shear deformation. We show that after several loading-unloading cycles 
the system reaches a non-equilibrium steady state (plastic shakedown) where it 
exhibits critical behavior characterized by the power law statistics of avalanches. 
The associated spatial fractality is revealed through the study of spatial correla- 
tions. A study of the power spectrum of the energy fluctuations shows the char- 
acteristic 1// noise which is another signature of scale free behavior. Overall, 
the critical manifestations of the plastic flow revealed by our simplified model- 
ing are in full agreement with experimental observations and with previous more 
elaborate numerical studies in 3D. 

In an attempt to understand the origin of criticality in analytical terms we per- 
formed a reduction of the system of ordinary differential equations describing our 
2D visco-elastic lattice model to a threshold automaton with conservative dynam- 
ics. As in ID case, we minimized out elastic fields and obtained an automaton 
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model for a discrete variable which serves as an indicator of the number of dislo- 
cations passing a given point. Despite partial linearization implied by the use of 
piece-wise quadratic potential and the replacement of time-continuous dynamics 
by a sequence of jumps, the automaton model exhibits the same critical behav- 
ior as the original ODE model. Our numerical study suggests that the discrete 
automaton behind the FK model has a statistical (weak) Abelian property which 
makes it amenable in principle to rigorous mathematical treatment despite the 
nonlocal nature of the discharge rules. 

A short announcement of some of the results of this paper can be found in 
[87]. 

3. The model 

Classical phenomenological theory of plasticity is applicable in both crystal 
and amorphous settings because it circumvents an explicit reference to the nature 
of the defects carrying inelastic deformation. Instead, it operates with a concept 
of an averaged non affine deformation represented locally by a tensorial measure 
of plastic strain. The evolution of the plastic strain is governed by a dissipative 
potential which is assumed to be a homogeneous function of degree one making 
the dissipation rate independent. Other internal variables characterizing harden- 
ing, softening, additional inelastic rearrangements, etc. are often added as well, 
each equipped with its own dynamics. For instance, in crystal plasticity the effec- 
tive plastic strain rate is parameterized by slip measures associated with different 
slip directions and evolved by coupled rate independent mechanisms. Likewise, in 
transformational plasticity of shape memory alloys the evolution of the phase/twin 
fractions is usually assumed to be governed by (coupled) rate independent dy- 
namic equations. 

Despite their success in fitting the observed memory behavior, the phenomeno- 
logical models of plasticity do not describe adequately the internal dynamics re- 
vealed by endogenous fluctuations because they do not deal explicitly with micro- 
scopic time and length scales. In particular they miss temporal intermittency and 
spatial fractality. 

An alternative program of describing plasticity directly in terms of defects re- 
sponsible for the flow is far from being complete despite many important recent 
successes [123, 55, 124, 67, 122, 121, 120]. First of all, in order to handle micro- 
scopic dynamics of dislocations one needs to go back to the elastic framework. 
Since the strain field of a dislocation decays rapidly away from the core region, 
the corresponding long-range distortions can be adequately described within the 
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framework of linear elasticity. Therefore dislocations are modeled as (line) sin- 
gularities of linear displacement fields bounding surfaces that carry the slip. The 
latter is represented in the elastic models as a quantized displacement discontinu- 
ity and this is the only place where the inherent discreteness enters the picture. 

The localization of dislocations in space and time requires constitutive as- 
sumptions which are outside the realm of linear elasticity. Such assumptions are 
usually presented in the form of additional kinetic relations of phenomenologi- 
cal nature describing the associated mesoscopic, rate dependent dissipation [125]. 
An additional problem in this approach is that independent phenomenological as- 
sumptions are also required to describe nucleation, annihilation and other topo- 
logical transitions. Notwithstanding all this uncertainty, the simulation of realistic 
3D flows requires tracking of a huge number of interacting defects which remains 
prohibitively expensive in terms of computational time. 

In search for an alternative description we first recall that a simple hybrid 
discrete-continuous dislocational model can be formulated by using the Peierls- 
Nabarro framework [10, 126, 106]. Suppose that the bulk elastic energy is defined 



where u{x) is the displacement field and the function is quadratic. The field 
u{x) is allowed to have discontinuities [u] which are penalized by surface energy 
of the form 



In Peierls-Nabarro model the function is assumed to be periodic and the evo- 
lution of the field u is governed by an appropriate gradient flow dynamics. The 
main technical difficulty in this approach is to track the a priori unknown surfaces 
of displacement discontinuity S{t) with arbitrary complex geometry. 

In order to avoid the explicit 'tracking' of the slip surfaces one can consider 
their 'capturing' by incorporating discreteness already into the bulk energy term. 
To this end one can introduce an internal length scale a and consider a mesoscopic 
energy which we schematically represent as 



Now u is the lattice field and the function / is periodic with infinite symmetry 
group [128, 111] . The discrete energy / contains information about both the bulk 
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energy 0(Vm) and the surface energy the mechanism of such splitting in 

the continuum limit was discussed in the context of fracture in [129] and in the 
context of plasticity in [127]. It is important to keep in mind that the discrete 
energy <l>(n) describes multi-particle elastic units and is therefore different from 
the conventional molecular dynamics potentials dealing directly with interatomic 
interactions. 

After the energy density in (1) is specified, one can formulate the rate depen- 
dent dynamics for the elastic fields carrying dislocations. Under the assumption 
that such dynamics is viscous and overdamped we obtain the following system of 
ODEs 

vu = -d^{u)/du, (2) 

where v is the ratio of the internal time scale and the time scale of the driving 
[109, 8]. Notice that instead of classical visco-elasticity we consider here an 'en- 
vironmental' viscosity which allows us to avoid certain dynamic degeneracies in 
the system without quenched disorder [101]. Notice also that the system (2) is not 
autonomous because time enters through the equations for boundary units as, for 
instance, in the case of cyclic driving in a hard device. 

The numerical solution of (2) may be rather challenging under generic load- 
ing. However, in the limit of quasistatic driving z/ — )■ 0, which can also be in- 
terpreted as the limit of zero viscosity, the dynamic system (2) can be replaced 
almost everywhere (in time) by the equilibrium system [8] 

d^{u)/du = 0. (3) 

The dynamics is then projected onto the branches of metastable equilibria which 
are defined as sequences of local minima (of the elastic energy) parameterized by 
the loading parameter. The crucial observation is that similar to the case of relax- 
ation oscillations, the dissipation-free phases of such dynamics are necessarily in- 
terrupted by the branch switching events which are instantaneous at the time scale 
of the quasi-static driving. Such reduction of an original dynamic problem opens 
the way towards minimizing out the elastic variables and reformulating the prob- 
lem in terms of a discrete integer valued variable labeling individual metastable 
branches [109, 8, 101]. 

In the next section we illustrate this general approach by using the simplest one 
dimensional setting where our goal will be modest: to make the idea of marginal 
stability fully transparent and to obtain the simplest examples of fluctuating plastic 
flow. Then we formally expand the model to the second dimension and show 
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that the ensuing minimal geometrical complexity is already sufficient to generate 
realistic fluctuational behavior. 

4. One dimensional setting 

The very first model of plasticity accounting for fluctuations was probably 
Prandtl's zero-dimensional model describing a particle driven through a spring in 
a periodic landscape [6]. This model involves marginal stability, captures hys- 
teresis, and reproduces rate independent dissipation [8] but misses the crucial 
interactions among the defects. To obtain criticality in this model one needs to 
considerably complicate the effective landscape, replacing, for instance, a peri- 
odic potential by a Brownian motion type potential [84]. Such fine 'tuning' of 
the interaction potential places this model in the class of phenomenological mod- 
els representing a natural development of the classical rheological models used in 
engineering plasticity. 

The simplest nontrivial 'first principle type' model with the energy (1) de- 
scribes the transformational plasticity of shape memory alloys [130, 131, 132, 
133, 8]. In this model each elastic element interacts with its nearest neighbors 
(NN) only through the total stress which then plays the role of a mean field. 
Having the advantage of utmost simplicity, this model possesses a nonphysical 
permutational invariance which introduces undesirable degeneracy. Therefore we 
subsequently augment this model by introducing linear interactions between next 
to nearest neighbors (NNN) and allowing them to be either attractive and repul- 
sive [102, 103]. Finally, we replace a bi-stable NN potential by a periodic one and 
introduce quenched disorder. 

4.1. NN model 

Consider a chain of particles forming in the reference configuration a regular 
ID lattice. Denote the reference particle positions by x° = z where we imply 
that the atomic spacing is taken as the length scale (a = 1). Introduce discrete 
displacement field Ui{t) and the discrete strain field ei{t) = Ui+i{t) — Ui{t). Our 
Fig. 1 shows the simplest representation of this mechanical system. In plasticity 
framework one should rather view the chain as a deck of rigid cards interacting 
through nonlinear shear springs and view the displacement field Ui{t) as being 
perpendicular to the direction of the coordinate axis. 

For the NN model the equations of motion (2) take the form 

viii + a(ui+i - Ui) - a{ui - Ui-i) = 0, (4) 
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Figure 1: A schematic representation of a one dimensional NN chain. 
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Figure 2: Stress-strain (a) and energy-strain (b) relations in the NN model with = 32, k = 1, 
V = 10"''. In (b) dashed line represents loading, solid line-unloading. 

where cr(e) = df /de is the elastic force resulting from stretching of a spring. We 
assume that the system is placed in a hard device with 



where t is the dimensionless time playing the role of the loading parameter. 

One can see that the total energy in the NN model is the sum of the energies of 
individual springs. This means that the springs interact only through the constraint 
imposed by the hard device. Indeed, in equilibrium cr(ej) = a, where a is the 
stress common to all springs; such infinite range of interactions is characteristic 
of paramagnetism. 

To facilitate the subsequent elimination of elastic variables we assume that the 
double well potential describing individual springs is piece-wise quadratic: 



Here k is the elastic modulus which we assume to be the same in both phases. The 



Uo(t) = 0,nAf(t) = t 




K 
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interval ej < 0.25 will corresponds to Phase 1 with equilibrium strain e = and 
the interval > 0.25, to Phase 2 with equilibrium strain e = 0.5. 

Our Fig. 2(a) shows the macroscopic strain-stress curve obtained from the 
numerical solution of (4) with v sufficiently small. As we load the system from 
a homogeneous state in Phase 1, the phase transition starts at the critical strain 
e = 0.25, where only one spring changes phase. The transformation in one spring 
corresponds to a stress drop as the other springs relax. Each subsequent stress drop 
is also associated with only one spring changing phase during each jump. One can 
see that the resulting macroscopic stress strain curve forms a wiggly plateau with a 
rudimentary intermittency manifesting itself through equally spaced oscillations. 
This behavior was studied in [133, 132] where it was shown that the individual 
springs can change phase only consequently. 
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Figure 3: Blow up of a schematic strain-stress relation for the NN model in the limit — > 
{N = 8). Solid lines show the path followed by the system during quasi-static loading. Thin 
lines are metastable branches. The black dots mark marginally stable states while the white dots 
correspond to local minima where the system stabilizes after a jump. 

As we mentioned in the Introduction, our model can be simplified further in 
the case of quasi-static loading which for the the dynamical system (4) corre- 
sponds to a zero viscosity limit . We begin with rewriting the piece- wise quadratic 
potential (5) as 

/(e.)= 1{^^-^^f■ (6) 

Here, we introduced a double-valued spin variable nii describing the phase state 
of a given spring with = for < 0.25 and nii = 0.5 for e.j > 0.25. At a 
given spin field, the displacement field Ui must satisfy the following set of linear 
equations 

{ui+i + - 2ui) = {rrii - rui^i), (7) 
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where the right hand side plays the role of 'dislocation density' [76]. If the phase 
field is known, the tri-diagonal matrix on the left may be inverted and the strain 
field can be obtained explicitly. Moreover, due to permutational invariance of the 
problem, one can always assume that there is only one defect and therefore the 
equilibrium branches can be parameterized by a single variable p = {1/N) ^ rrii 
which changes in the interval < p < 0.5 and describes the 'phase fraction' in 
the mixture (see [132]). 

Several metastable branches parameterized by p are shown schematically in 
Fig. 3. The first branch with all springs in the first energy well (Phase 1) cor- 
responds to j9 = 0. It ends at point a where this homogeneous state becomes 
unstable. Suppose that we drive the system by quasi-statically changing the total 
strain. Then at point a one spring changes phase and the system undergoes a dy- 
namic transition to the new local minimum identified in Fig. 3 by point (3. This 
fast process is accompanied by a jump in stress; as a consequence of this jump a 
finite amount of energy dissipates into heat. 

From point (3 the driven chain evolves along the second equilibrium branch 
corresponding to p = 0.5 /N which in turn becomes unstable at point 7. Here 
again one spring switches phase, the stress drops and the system finds itself in a 
new local minimum identified as point 6 on the branch with p = 1/N. The process 
continues till all springs change phase and the system reaches the second homo- 
geneous branch corresponding to p = 0.5. Interestingly, the way we introduced 
viscosity ensures that the transformation proceeds as a single front which is not 
the case if we directly deal with the degenerate problem at = or if we consider 
classical visco-elasticity [101]. 

A full analysis of the thermodynamics of the NN model including the compu- 
tation of the 'heat to work ratio' can be found in [8]. 

4.2. Automaton model 

After the elastic fields are minimized out (either by the inversion of the matrix 
in (7) or by the Fourier transform as we show below), the ensuing problem in 
terms of the variables rrii can be reformulated as a spin-valued automaton. 

To be more precise suppose again that initially all the units are in the Phase 
1, meaning that the phase variable vanishes everywhere. Since the field rrii is 
known, one can compute Cj by solving (7). Then one can find a and locate the 
system on the equilibrium branch with p = 0. The stability condition in this case 
reduces to an inequality a{p, t) < cr'^, where a'' = 0.25 is the upper spinodal stress 
[8]. Suppose that the stability condition is violated at time ti. Then we need to 
update the phase configuration p = p + 0.5 /N which leads to a stress drop from 
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a{p,ti) to a{p,ti). Notice that the choice of the flipping element nii remains 
arbitrary in the automaton why in the ODE framework it is uniquely prescribed 
by the dynamics. To overcome this (artificial) non-uniqueness in the automaton 
setting we later consider two regularizing mechanisms: NNN interactions and/or 
quenched disorder. 

Along the new metastable branch p = p one can continue computing cr(p, t) 
till the spinodal stress is reached again and then it is necessary to make the next 
update. The stress oscillations continue till we reach the homogeneous Phase 2 
with p = 0.5. During unloading the stability condition changes into d-(p, t) > ac, 
where ac = is the lower spinodal stress and the update rule becomes p — )■ 
p — 0.5/N. The homogeneous branch corresponding to Phase 1 is reached again 
when p = 0. We do not present here the strain-stress curve obtained from the 
automaton-based simulations because it is practically indistinguishable from the 
results obtained by solving directly the dynamic equations (4) with v = 10~^ 
shown in Fig. 2. 

In order to show that on the yielding plateau the system evolves through a se- 
quence of marginally stable states, we plot in Fig. 2(b) the energy of the system 
against the applied strain for one loading-unloading cycle. When the system is in 
one of the pure phases (affine or Cauchy-Bom state) it resides in a global mini- 
mum of the energy until the Maxwell threshold aM = is crossed (during either 
loading or unloading). Then the global minimum becomes a local minimum and 
eventually the homogeneous state becomes unstable as the system reaches either 
upper or lower spinodal stress. After that the system evolves through the sequence 
of non-affine states where different springs occupy different phases. Notice, how- 
ever, that in this non-affine stage the stress in the chain never deviates considerably 
from one of the spinodal stresses (see points a, (3, 7, 5, etc. in Fig. 3). It is not 
hard to see that in the continuum/thermodynamic limit N ^ 00 the jumps become 
infinitesimal and the ensuing yielding plateau collapses the spinodal [8]. There- 
fore, in the continumm limit the yielding system is marginally stable on the whole 
yielding plateau. 

In summary, plasticity in this simplest model is associated with elastic states 
of minimal stability. This peculiar feature of plastic flow may be considered as a 
first manifestation of criticality. In addition, the model exhibits stress drops that 
can be interpreted as avalanches. The problem is that these avalanches are all of 
the same size and the corresponding microscopic transformation events are either 
un-correlated (zero viscosity limit) or over-correlated (finite viscosity case). 
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4.3. NNN model with 'ferromagnetic' interactions 

As we have seen, the NN model lacks the short range interactions capable 
of generating extended dislocation cores. To bring the missing internal length 
scale into the theory one can add to the NN model an additional 'ferromagnetic' 
interactions favoring local homogeneity of the strain field [134]. 

The simplest way to introduce short range interactions is to consider an NNN 
model with the energy density [102] 

f{ed= ^(e.-m,f + + («) 

Here nii is the spin variable introduced in the previous subsection and A is the 
strength of short range interactions (we consider only the automaton version of 
the model). The 'ferromagnetic' interactions ensuring the localization of phase 
boundaries corresponds to A < 0; the case of 'anti-ferromagnetic' interactions 
A > 0, favoring the formation of fine mixtures, will be considered in the next 
subsection (see also [135]). 

In the quasi-static limit discussed above, the equilibrium equations for the 
NNN model can be written in the form 

[ui+i + Ui-i - 2Mi] + (A/ k) [ui+2 + Ui-2 - '2ui] =mi- rrii^i, (9) 

where again the 'dislocation density' in the right hand side plays the role of a 
source of elastic strain. Due to non locality of the NNN model, we need to com- 
plement the hard device boundary conditions used in the NN model by two ad- 
ditional conditions. One way to avoid surface boundary layers is to add two fic- 
titious springs just outside the boundaries of the chain ensuring the periodicity 
conditions e(0) = e(l) and e(A^ + 1) = e(A^). This leads to additional boundary 
terms in the energy, however, the elastic problem can be again solved explicitly 
[102]. The remaining problem for the spin field mj(t) can be formulated as an 
automaton which is similar to the one described in the previous section and we 
omit the details. 

The macroscopic strain-stress curves for the NNN model with A < are 
shown in Fig.4(a). In addition to the horizontal plateaus, that again correspond 
to the regimes when an isolated phase boundary sweeps through the chain (propa- 
gation), one can also see two peaks at the onset of the direct and inverse transfor- 
mation (nucleation) and two smaller anti-peak describing the disappearance the 
existing phase boundaries (annihilation). In the computational experiment shown 
in Fig.4(a) four springs out of 24 participate in the nucleation event while the 
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Figure 4: (a) Strain-stress relation for the NNN automaton, (b) Evolution of the total energy. Here 
and in other similar plots N = 32, k = 1 and A = —0.01 

propagation involves successive jumps of individual springs. The necessity of the 
peak in the NNN model and the nature of the nucleation event in the continuum 
limit were discussed in [102]. 

In Fig. 5(b)) we show the macroscopic energy-strain relation for this model 
(loading only). In contrast to what we have seen in the NN model, here the en- 
ergy increases in the systematic manner as the system evolves along the plastic 
plateau. This means that plastic dissipation is accompanied by an additional en- 
ergy build up which can be vaguely qualified as the 'cold work' [8]. The origin of 
this effect in our model is the particular form of the NNN term in the energy (8) 
which, in addition to penalizing gradients, also contributes to homogeneous elas- 
ticity. Interestingly, similar continuing growth of the energy along the flat-stress 
yielding plateau have been observed in several microscopic (MD-type) models of 
amorphous plasticity [172]. 

In order to avoid such hidden contributions to bulk elasticity, one can intro- 
duce short range interactions in a different form. The simplest way to keep the 
ferromagnetic coupling intact while removing the NNN contribution to homoge- 
neous elasticity is to introduce a three-particle interaction mimicking the strain 
gradient (SG) term in the discrete setting 

m= ^(e,-m,)^ + ^(e,-e,_0'. (10) 

It is not hard to see that in the case (3 > this energy is equivalent to the NNN 
energy (8) with A < up to a 'local' term depending on only. This means that 
the 'cold work' is eliminated in the SG formulation which is confirmed by our Fig. 
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Figure 6: (a) Strain-stress relations for the NNN model with A = 0.01, other parameters are the 
same as in the other similar graphs. A step in the center of the plateau corresponds to elastic 
deformation of a 'binary' phase with lattice scale oscillations; (b) Evolution of the total energy. 



5(b). Notice that outside the nucleation (annihilation) events the behavior of the 
NN and SG models are quite similar with isolated phase boundaries propagating 
along the chain as the total strain increases. The main difference is that in the NN 
model a phase boundary (mimicking dislocation) is atomically sharp, while in the 
SG model its core has a finite size. 

Despite all these new features, neither NNN nor SG model are successful in 
predicting the realistic statistical structure of avalanches. Indeed, both models 
exhibit one correlated 'snap' event (nucleation) and a succession of equally sized 
'pop' events (propagation). It is clear that this fluctuation/avalanche structure is 
still too simple to be compared with realistic experiments. 
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4.4. NNN model with 'anti-ferromagnetic' interactions 

As it is well known, the nonlocal kernel describing linear elastic interactions 
in dimensions 2 and 3 has both ferromagnetic and anti-ferromagnetic contribu- 
tions [136, 137, 138]. Therefore it is of interest to see how plastic behavior in 
our toy model changes if 'ferromagnetic' NNN interactions with A < are re- 
placed by 'anti-ferromagnetic' NNN interactions with A > 0. In the ID setting, 
anti-ferromagnetic interactions favoring formation of small scale microstructures, 
should be understood as a poor man's attempt to capture the gradient nature of 
the strain (elastic compatibility) and the effect of the constraints placed by higher 
dimensional boundary conditions [137]. 

Our numerical results for the anti-ferromagnetic NNN model driven in a hard 
device (automaton version) are summarized in Fig. 6. The most important feature 
is that in addition to two homogeneous phases the model exhibits the third phase 
representing binary lattice scale mixture which is homogeneous only in average 
(see Fig. 7). This new phase shows binary oscillations with neighboring springs 
occupying different phases; the energetic preference of such non Cauchy Bom 
(non-affine) phase arrangements follows from the studies of the global minimum 
of the elastic energy, e.g. [139, 140, 141, 103, 142, 135]. 




X X 



(a) (b) 

Figure 7: Strain profile evolution in the driven NNN model with A = 0.01. (a) Propagation of the 
binary phase along the first half of the yielding plateau; (b) homogeneous binary phase with lattice 
scale oscillations on the elastic step separating two yielding regimes in Fig. 6. 

The salient feature of the macroscopic response in the anti-ferromagnetic NNN 
model is the disappearance of the nucleation peak. One can see that on the first 
half-plateau of the strain-stress curve (Fig. 6(a)), the homogeneous phase is pro- 
gressively replaced by the growing 'binary phase' and the corresponding front is 
shown in Fig. 7(a). Once again, each drop in stress profile corresponds to one 
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Figure 8: (a) Strain-stress relation in the driven SG model with /3 — — 0.2, (b) evolution of the 
total energy at varying applied strain. 

spring changing phase. Notice that the thresholds for the nucleation of the binary 
phase and for its propagation are similar. By the end of the first half -plateau the 
binary phase occupies the whole domain and the step in the middle of the strain- 
stress curve (Fig. 7(b)) describes purely elastic deformation of the binary mixture. 
Along the second half -plateau the mixture phase is gradually replaced by the sec- 
ond phase and we again observe a propagating front. The dynamics of such fronts 
has been studied in [143, 144]. 

From Fig. 6(b) we see that 'ferromagnetic' model, the bulk contribution to 
the energy in the 'anti-ferromagnetic' model changes its sign. The main plastic 
phenomenology, however, remains the same even though the microscopic spatial 
arrangement of phases changes from phase separation to fine scale mixing. 

For completeness, we also present here the results of the simulations for the 
anti-ferromagnetic SG model with /3 < 0, see Fig. 8. The strain profiles in 
this model are basically the same as in the NNN model with anti-ferromagnetic 
interactions (A < 0), however, the energy profile shows a new feature: a negative 
slope in the first segment and a positive slope in the second. This implies that the 
negative contribution to the cold work along the first half-plateau is compensated 
by the equal in magnitude positive contribution on the second half-plateau. 

To summarize, the ID anti-ferromagnetic NNN model and its SG analog de- 
scribe rate independent dissipation and show a new feature which can be inter- 
preted as a propagation of regular slip band microstructures. Still, both models 
predict highly coherent fluctuations structure and therefore none of them captures 
the complexity of the experimentally observed spatial and temporal correlations. 
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4.5. Quenched disorder 

In modeling of realistic material response one cannot neglect disorder due to 
various chemical and mechanical imperfections at the level of a lattice. Such dis- 
order can influence the macroscopic behavior of the system by facilitating inho- 
mogeneous nucleation of defects on one side and by creating a variety of pinning 
obstacles for the propagating defects and defect microstructures, on the other side. 
In this subsection we show how the quenched disorder influences the macroscopic 
response in the basic NNN model with ferromagnetic interactions. 

A disorder can be incorporated into the energy (10) in different ways, for in- 
stance, by randomizing the size of the barriers separating individual energy wells 
[133]. In this paper, due to technical reasons, we follow a different approach 
originating from the RFIM model [145] where disorder is represented by a ran- 
dom field biasing different phases in different spatial points. The ensuing RSSM 
model (random soft/snap spring model [1 10]) is characterized by the energy 

/(^i) = - ^if + + ^i-if - hiei, (11) 
where A < and are independent random numbers distributed according to 

P{x) = {V2^y^exp{-{xf/{2a^)). 

In all our numerical illustrations we put cr^ = 0.01; we have checked that the 
value of dispersion does not affect our statistical results even though it affects 
some macroscopic features, for instance, the degree of hardening (see [101]). 
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Figure 9: Strain-stress relation for the NNN ferromagnetic model with quenched disorder. Here 
~ 2048, other parameters are the same is in similar graphs shown before. 
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A typical strain-stress relation for a particular realization of disorder is shown 
in Fig. 9. We observe that the nucleation peak persists, however the successive 
stress jumps due to phase boundary advances have random amplitudes and each 
avalanche now involves different number of springs. Between the jumps the sys- 
tem is trapped/pinned in metastable configurations where it evolves purely elas- 
tically. The instabilities can be again associated with a succession of minimally 
stable states where at least one element have reached the spinodal state. Notice 
that the disorder in this model does not evolve and therefore the shakedown state 
is reached already during the first cycle. 




(a) (b) 



Figure 10: (a) Time series for the dissipated energy E{t) during one loading cycle in the NNN 
model with quenched disorder shown in Fig. 9 ; (b) blow up of three successive avalanches. 

We observe that as in the case without disorder, each jump at time t is asso- 
ciated with a finite dissipation which can be computed as the energy difference 
before and after the jump 

E{t) = Y,iMt+)-Mt-))- (12) 

i 

At a given t the value of E(t) may be either zero (elastic deformation) or positive 
(irreversible jump). In a system with disorder a single jump, originating from a 
phase change in one unstable element, does not necessary lead to a stable state and 
therefore an avalanche may contain several elementary jumps. We have already 
encountered this phenomenon while dealing with nucleation requiring simultane- 
ous transition in several elastic elements. In the case of nucleation, however, we 
had to deal with massive transformation at one value of the loading parameter. In 
a disordered system typically only one element becomes unstable at a given value 
of the load, however, several such jumps may take place before the next elastic 
stage is reached. 
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A typical time series for the dissipated energy E{t) is shown in Fig. 10(a). 
We observe that jumps are separated by quiescent time intervals where E(t) = 0. 
The shape of several typical avalanches is shown in Fig. 10(b) where we do not 
see any significant variation of time and length scales. Notice also that avalanches 
are very short, corresponding to 10-20 time steps. All this indicates that despite 
marginal stability of the system we record a succession of local events which fail 
to trigger global rearrangements. 
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Figure 1 1 : Log-log plot of cumulative probability function (CPF) for the dissipated energy E in 
the NNN model shown in Fig. 9 and Fig. 4.5 . Dashed line is a gaussian fit: P{x) = Aexp{—{x — 
/^)V(2cr2)), ^jjjj ^ ^ 8.9078, fi = -0.0071 and a = 0.0034. 

The total dissipated energy during an avalanche separated by two silent inter- 
vals can be written as 

A 

where summation is over the avalanche duration. The cumulative probability dis- 
tribution [146] of the variable E is shown in Fig.ll. One can see that the dis- 
tribution is close to Gaussian if we exclude few large scale events corresponding 
to macroscopic nucleation/annihilation. The fact that we do not record a power 
law signal implies that this ID system does not exhibit criticality and that outside 
highly synchronized nucleation events the avalanches remain largely uncorrelated. 

4.6. Explicit elimination of elasticity 

In the case when the elastic potential is piece-wise quadratic, one can solve the 
linear elastic problem analytically and present the automaton explicitly in terms 
of spin variables. Given our periodic boundary conditions, it is natural to invert 
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the linear elasticity operator in the Fourier space. The discrete Fourier transform 
of a lattice field Xk is defined by X{q) = l/Nj^k^ke''"'', where q is the ID 
wave vector with values q = 2Txj/N , j = 1, ...N. The inverse Fourier transform 
is given by Xk = Yuq^i^Y"^^- 

Consider first the NN model. By transforming (7) with added quenched disor- 
der into Fourier space we obtain 

Kc{q)u{q, t) - s-{q){Km{q, t) + h{q))) = 0, (13) 

where 

c{q) = 2[cos(g) - 1], 
SxiQ) = (1 - cos(g) +isin(g)). 

The field rh{q,t) is the Fourier transform of the discrete spin field rrii and the 
field h{q) is the transform of the random field hi \ both fields represent sources of 
elastic strains. The difference between these two sources is that the field h{q) is 
fixed (quenched disorder), while the field m(g, t) is evolving in accordance with 
the automaton rules which we formulate below. 

First, we decompose our linear elastic field into a homogeneous part associ- 
ated with prescribed time dependent affine deformation on the boundary and an 
inhomogeneous part which is a solution of a problem with periodic boundary con- 
ditions. More precisely, we assume that the g = component of the strain field is 
controlled by the loading while the g 7^ component can be obtained by solving 
a periodic problem with given sources. 

The expression for the non affine component of the displacement field can be 
obtained from (13) 

^^^^^^^s-Aq){nrn{qt)^h{<l))_ 
Kc[q) 

The corresponding strain field e^^{q) = s^{q)u{q) where 

= -(1 - cos(g) - i sin(g)), 

can be written as 

6^^(g) = ^i^^^(«:m(g) + h{q)). (15) 
Kc{q) 

The affine component of the deformation field eo = J2i fully defined by 
the loading conditions 

eo(g) = 5iq)t. (16) 
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The total strain is then equal to 



e(g) = eo(g) + e^^(g). 



The same method can be applied to the NNN model and the SG models with the 
replacement of e^^{q) by either e^^^(g) or e'^'^(g). We obtain 



and 



Kc{q) — X sin^(g) 



[Km{q) + h{q)). 



Km{q) + h{q)). 



(17) 



(18) 



Kc{q)-mQW 

Such elimination of the elastic fields allows one to formulate the condensed au- 
tomaton model for the spin variable rh(q) in explicit form. 

4. 7. SG model with periodic potential 

As an application of the reduction method presented in the previous section 
we consider here a discrete chain with periodic NN potential. Once again we 
circumvent dynamic simulations and move directly to the study of the automaton 
model corresponding to the inviscid limit z/ — )■ 0. 




Figure 12: Periodic potential and its piece-wise quadratic approximation. One elastic domain is 
shadowed. 

We chose the periodic potential to be piece-wise quadratic (see Fig. 12) and 
assume that in one period 



(19) 
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where d = 0, ±1, ±2, ... is a new integer-valued spin variable which describes 
quantized slip; each value of d defines a quadratic energy well with a period (2^°). 
At di field given, the elastic strains Cj must satisfy 

{d, - l/2)(2e°) < e, < {d, + l/2)(2e°). (20) 

To regularize the NN model we add to the elastic energy the ferromagnetic SG 
term as in (10) and introduce a quenched disorder as in (1 1). The elastic field can 
now be decomposed into three components 

e, = eo + ef + ef = {6o(g)}r' + {e\q)};' + {e\q)};'. (21) 

Here the affine component eo{q) is given by (16); the Fourier transform of the 
residual strain due to quenched disorder can be written as 

e\q) = L{q)h{q). (22) 

where the Green's function is 

and finally the Fourier representation of the strain field due to plastic slip is 

e\q) = L{q)d{q). (24) 

Now, according to (20), a metastable configuration with a given distribution di is 
defined (is stable) within the following limits 

-e-t- {e\q)}7' < G,{d) <e-t- {e\q)}-\ (25) 

where 

did) = {L{q)diq)}-' - i2e)d^. (26) 

When Gi{d) reaches one of the thresholds, the integer parameter di must be up- 
dated 

d^^di + Mi{d), (27) 

where 

+ 1, ifG,{d)>e-t-{e\q)}-\ 
M,{d) = <; - 1, if G.(rf) <-e-t- {e\q)}-\ (28) 
otherwise. 
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Figure 14: Stress-Strain relations in two models with periodic potential and quenched disorder 
(a)/3=0 (NN model) and (b)/3 = 0.01 (SG ferromagnetic model). Inserts show the corresponding 
strain profiles. 



In a driven system an initially stable distribution of slip di can get destabilized 
because the thresholds evolve with time. The update of the slip di given by (28) 
leads to a redistribution of the elastic strains ef. In Fourier space the corresponding 
'discharge rules' take the form 

e(g)->e(g) + L(g)M(g). 

In Fig. 13(a) we show the physical image L{x) of the discrete kernel L{q) for 
/3 > (ferromagnetic SG model) and (3 < (anti-ferromagnetic SG model). One 
can see that if an elementary positive slip takes place in an element, the elastic 
strain is increased in several neighboring elements and this 'influence domain' 
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increases with absolute value of /3. In NN model where /3 = the strain of only 
the unstable unit increases while the strain of all another units decreases which 
is a typical mean field behavior. In the regularized model with /3 7^ the local 
response is more complex as one can see in Fig. 13(b). 

One of the crucial points in the automaton formulation is the choice of the 
update strategy when several elements become simultaneously unstable. To avoid 
this problem we introduce slip dependent disorder hi{d) which signifies that ran- 
dom bias fields may be different in different energy wells. Physically this means 
that imperfections may affect unevenly the repeated slips along the same slip di- 
rection. Similar assumption is usually made in the meso-scopic models of fric- 
tion/amorphous plasticity where new random thresholds are assigned to the newly 
formed bonds, e.g. [28]. 

In Figs. 14(a) and 14(b) we show the strain-stress relations in the periodic 
problem with (3 > (SG ferromagnetic model) and (3=0 (NN model). In both 
cases, the macroscopic response is very similar: deformation is irreversible, springs 
always remain close to the marginal stability limit and the transformation proceeds 
intermittently through random avalanches. However, the strain profiles ej in the 
two models are different as shown in the inserts in Figs. 14(a) and 14(b). Thus, in 
the local model with (3 = 0, the strain field is chaotic at the lattice scale due to the 
presence of disorder. Instead, in the nonlocal model which penalizes interfaces we 
see a more orderly (but still random) arrangement of the domains with different 
level of slip. In both cases the statistical distribution of dissipated energy during 
avalanches remains Gaussian as in the case of transformational plasticity (a model 
with a double well potential). 
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Figure 15: Power spectrum in the ID automaton model with periodic NN potential and SG ferro- 
magnetic interactions. 
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It is also of interest to analyze the power spectrum of the signal E(t) [45, 147] 



where T is the total duration of the time series, / is the frequency. The log-log 
plot of the power spectrum for the SG model with /3 = 0.01 is shown in Fig. 4.7. 
We see the signature of a typical random process when the energy distribution is 
constant over a broad range in the low frequency domain. The rapidly decaying 
part of the spectrum at high frequencies is just a sign that the exerted power is 
finite. This behavior suggests a short range of temporal correlations which is 
incompatible with criticality. We recorded similar structure of the power spectrum 
in all other ID models, with and without short range interactions. 

5. Two dimensional setting 

Even though our driven ID models exhibited marginal stability, we could not 
detect in our numerical experiments any signs of scale free correlations. This may 
be associated with the simplicity of the geometrical structure of interaction which 
prevents local events from triggering global rearrangements. As we have seen, 
this conclusion remain valid even in the presence of strong Gaussian disorder that 
is not sufficient in the ID setting to generate avalanches with power law statistics. 
It is then natural to try to augment the model by introducing an additional spatial 
dimension. 

5.1. The model 

The simplest 2D extension of the models presented in the previous sections 
can be viewed as a series of ID NN chains linked trans versally by linear elastic 
springs (see Fig. 16 ); the ensuing 2D lattice model can also be interpreted as an 
array of coupled FK chains [148, 149]. Notice that in 2D the short range NNN 
or SG terms are not necessary because the nonlocal elastic interactions take place 
already in the NN model by virtue of an additional geometrical dimension. We 
also emphasize that while in ID we could discriminate only homogeneous slips 
(see Fig. 17), in 2D one can also describe the structure of the dislocation cores (see 



We assign to each particle with coordinates (z, j) defined on a x square 
lattice (a = 1) a horizontal scalar displacement field Uij. We associate with linear 
longitudinal springs the strain measure 6', ^ = Mi+i,fc — Ui^k and with nonlinear 
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Fig.18). 
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Figure 16: The schematic representation of the 2D model. Shear springs AB and CD are described 
by periodic potentials, longitudinal linear springs AC and BD are described by quadratic potential. 
Displacement field is horizontal. 



shear springs - the strain measure ^ = Ui^k+i — Ui^k- The total energy of this 
highly anisotropic lattice can be written as 

Hu) = J2fie.,,^,,), (30) 

where the potential / is assumed to be quadratic in 9 and periodic in ^. More 
specifically, we set 

f{9,„ 6,,) = 9ii,j) + ^iO,,Y - hl^.j - hy,,. (31) 

where g is a periodic function (see Fig. 12), in particular, in our dynamic (ODE 
based) model we use 

^?(0 = (27r)-2(l-cos(27rO). (32) 

The two lattice functions hl'j are independent Gaussian random variables intro- 
ducing quenched disorder. 

Observe that in the limit K — oo, we must have = and our 2D model 
reduces to a one dimensional model with periodic potential which we have already 
discussed in the previous sections (see Fig. 17). When K < oo the variables 9 and 
^ are no longer independent and the corresponding long-range interactions can be 
revealed by minimizing out a linear variable 9 which is lacking in the ID model 
[136, 150]. 

An elementary slip in 2D can be illustrated already in the constrained model 
with = oo if the adjacent atomic planes are displaced by one or several lattice 
spacings (Burgers vector) without changing the energy, see Fig. 17. Within the 
classical interpretation, one usually relabels the nearest atoms after the slip takes 
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(a) (b) 

Figure 17: A quasi-one dimensional slip deformation in a 2D model: (a) with relabeling, (b) 
without relabeling; b is the Burgers vector 

place: notice that the bond in Fig. 18(a) between the atoms A and B is broken, 
therefore, a new bond has to be added between the atoms B and C, which are 
now nearest neighbors. This approach gives rise to discontinuous displacement 
fields and requires additional phenomenological hypotheses to deal with evolving 
discontinuities. In our model we do not perform the relabeling and consider com- 
patible representation of the lattice before and after slip, see Fig. 17(b). Notice 
that in this case, one cannot use linear elasticity and has to introduce a periodic 
potential which ensures lattice invariance with respect to finite shears. 

If K < oo the plastic flow does not occur by a slip of complete atomic planes. 
Instead the slipped atomic planes end inside the crystal forming dislocations. In 
Figs. 18(a) and 18(b) we illustrate the role of relabeling in the classical repre- 
sentation of dislocations and compare it with the compatible description of large 
shears adopted in this paper. 




(a) (b) 



Figure 18: A single dislocation in a 2D model: (a) classical model with relabeling or (b) compati- 
ble description without relabeling. 

To complete the formulation of the problem we assume that the lattice is sub- 
jected to a time dependent shear in a hard device. In our framework driving with 
constant strain rate can be presented as displacement controlled boundary condi- 
tion 

Ui^N-l = Uifi + t 



32 



for 2 = 0, — 1. This boundary condition can also be rewritten in terms of the 
overall shear strain as 

Af-l 

where t is again the slow time playing the role of loading parameter. In longitudi- 
nal direction we assume the periodic boundary conditions given by 

for j = 0,N — 1. The overall axial strain must then satisfy 

N-l 

J2 ^^^k = 0. 

i=0 

Observe that our elastic energy depends on two material parameters: K, rep- 
resenting the ratio of the bulk and shear moduli and the ratio of the Burgers 
parameter and the interatomic distance. In addition, the behavior of the system 
may be affected by the variance of the quenched disorder a. The natural question 
is where in the space of parameters one can expect to see the critical regimes. 
For large K the model becomes one dimensional, while for large the dislo- 
cations disappear. When K is small the model again becomes one-dimensional 
and for small we lose lattice trapping. At very large disorder, dislocation mo- 
bility is very limited (POP regime, [110]) and the critical behavior can hardly be 
expected. On the other hand, at small disorder, one can expect synchronization 
(SNAP regime, [110]). The detailed reconstruction of the corresponding bound- 
aries presents a formidable task which is outside the scope of this paper. A lim- 
ited parametric study shows robustness of the scale free behavior, in particular, 
we observed that the power law statistics was disorder insensitive in the interval 
0.01 - 0.2. 

5.2. The dynamic problem 

As we have seen in the ID case, a quasistatic evolution of the driven system 
can be studied either in a dynamic setting with small but finite viscosity or in an 
automaton setting corresponding to the limit of infinitely slow driving. In the 2D 
setting we first formulate the dynamical model with small viscosity and smooth 
periodic potential (32). Later we compare the results with the automaton model 
implying zero viscosity and operating with a piece-wise quadratic approximation 
of the periodic elastic potential. 
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To facilitate the Fourier representation we reintroduce the finite difference op- 
erators st and s,t 

Then the axial and shear stresses can be written as 

a- = Kstu,,, <J = 2n-' sin(27r(4u,,)), (34) 
and the equation of motion (2) takes the form 

lyiiij = Ks-al^^ + s-a'^l, (35) 

where 

^-n^^ - _ XX - xy _ xy xy 

In order to avoid numerically challenging non-local spatial derivative terms in 
(35), we use the Fourier method. The discrete Fourier transform of a 2D lattice 
function is defined by 



k I 

where q = {qx.Qy) = ^) is the wave vector. By mapping our dynamic 
equation (35) into Fourier space we obtain 

pii{ci) = Ks-{q)a^%q) + 5;(q)a^nq) - ^(q), (38) 
where we defined 

"Sa (q) = (1 - cos(ga) + ism{qa)), a = x,y, 
and added an inhomogeneity due to quenched disorder 

To solve (38) numerically we use a semi-explicit Euler time discretization. The 
axial strain a^^ is a linear term with respect to the displacement u and we can treat 
it implicitly in Fourier space. The non linear term a^^ is first calculated in real 
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space and then transformed into Fourier space and treated explicitly. As a result, 
we obtain 



where 



(q) = -(1 - cos(ga) - i sin(ga)), a = x,y, 

and At is the time step. After the field n*(q) is known, the strains can be computed 
in Fourier space by simple multiplication. To obtain the real space solution we 
need to invert the Fourier transform by using the formula 



^M = EE^(q)e'^'''^'"^- (40) 

k I 

Notice that the FFT method automatically enforces periodic boundary condi- 
tions for the displacement field Uij. For instance, the conditions J2iLo^ ^i,k = 
are automatically satisfied. Driving in shear J2fJo^ = t can be again im- 
plemented by means of decoupling of the homogeneous shear strain component 
eo(q) =t5(q). 

5.3. Numerical results in the dynamic problem 

In this subsection we present the results of our numerical experiments with 
implicit-explicit FFT method described above. In our tests we always take K = 2 
and = 512 while our initial state is always dislocation free. 



Figure 19: A pre-strained zone in the form of an ellipse representing a local inhomogeneity and 
triggering dislocation nucleation. 

As a first illustration we show in Fig. 19 the strain state due to internal pre- 
stress hjj = 0.1 inside an ellipsoidal inclusion. In the driven system such inho- 
mogeneously stressed regions play the role of nucleation centers for dislocational 



35 



Figure 20: A fragment of the deformed lattice with two dislocation dipoles nucleated under the 
applied shear strain around an imperfection shown in Fig. 19). Red and blue colors indicate stress 
concentration corresponding to dislocations of different signs. 

dipoles (mimicking dislocation loops in 2D). In Fig. 20 we show a fragment of 
the deformed lattice with two dislocation dipoles nucleated around an ellipsoidal 
imperfection as a result of the application of homogeneous shear strain on the 
boundary of the domain. The local shear stress given by dg{^)/d^, is indicated 
by colors: red, for positive shear stress and blue for negative shear stress. One 
can see that the stress reaches its maximum inside the four distinct dislocation 
cores. As the applied strain is increased the dislocations of different sign move in 
opposite directions. Eventually they reach the boundaries which means that two 
full atomic layer have slipped. 

In Fig. 21, we show a typical dislocational configuration in a crystal with 
random quenched disorder subjected to a finite shear. Once again the colors indi- 
cate the location of dislocations with different signs. Here we already see several 
atomic planes that have experienced either singular or multiple slip. Notice also 
that the dislocational dipoles have a tendency to concentrate around regions with 
strong inhomogeneity. The corresponding stress field at larger scale is shown in 
Fig. 22. We observe that dislocations, which interact through long-range elastic 
fields, organize themselves into complex microstructures which appear random 
but turn out to be highly correlated. In our simulations we see that in the steady 
state regime plastic activity reduces to intermittent dislocational exchanges be- 
tween stable clusters containing multiple dislocational dipoles. While these clus- 
ters remain mostly unchanged from one cycle to another, they have a finite lifetime 
as in experimental observations reported in [151]. 

In Fig. 23(a), we show the evolution of the dislocation density during the first 
six cycles of loading/unloading which we compute by localizing and identifying 
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Figure 2 1 : A fragment of a generic dislocational configuration during a steady state simulated plas- 
tic flow, red and blue colors indicate stress concentration around dislocational cores of opposite 
signs. 

individual dislocation cores. We observe an initial overshoot and the subsequent 
stabilization which indicates that the system has reached the shakedown state. In 
Fig. 23(b) we present the corresponding macroscopic strain-stress curves exhibit- 
ing plastic hysteresis. The size of the hysteresis loops decreases with cycling, 
which indicates that the rate of dissipation diminishes as the system get stabilized 
in the shakedown steady state. 

One can see that the macroscopic description of plasticity in our 2D model is 
quite realistic. At the microscale we see that the initial incipient inhomogeneity 
is augmented and modified during the cycling. This leads to the formation of 
statistically stable distribution of residual stress which is largely independent of 
the initial disorder. 

5.4. Statistics of avalanches in the dynamic model 

The fluctuations associated with intermittent dislocational dynamics are de- 
tected in experiment through the acoustic emission (AE) [21, 40]. The of intensity 
of the AE is accessed by measuring the square of the voltage which is expected 
to be proportional to the dissipated energy. Therefore by computing the rate of 
dissipation 




(41) 
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Figure 23; (a) Evolution of the dislocation density in the system subjected to cyclic loading. The 
inset shows small scale fluctuations; (b) Macroscopic strain-stress hysteresis during the first six 
cycles under strain controlled loading conditions. 
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Figure 24: (a) Time dependence of the dissipated energy during several cycles; (b) The same signal 
at shorter time scale (magnification of the highlighted region in (a)). 



in our simulations we can make comparison with experimental data. The typical 
time series for E(t) is shown in Fig. 24 at two different scales. One can see that 
the signal has a characteristic spiky appearance and that the complexity of the 
time series is not reduced with magnification which is an indication of scale free 
behavior. 

To separate individual avalanches we introduce an irrelevant threshold and 
define the avalanche energy by integrating the dissipation rate over the duration 
of the avalanche T: 

rt+T 

E = vN~^ / Kjdt- (42) 

The complexity of the intermittent signal is characterized by the probability dis- 
tribution of avalanches P{E) which in our case stabilizes after few first cycles. In 
Fig. 25(a) we show the stationary distribution of avalanches exhibiting a robust 
power law pattern for over 4 decades with exponent e ~ 1.6 ±0.05 obtained by the 
maximum likelihood method [146]. This value of the exponent is in perfect agree- 
ment with experiments in ice crystals and fits the generally accepted range 1 .4-1.6 
[35, 16, 42, 36, 35, 55]. It is also consistent with the value obtained for 2D col- 
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loidal crystals [152] and with the value obtained in phase field crystal simulations 
[68]. Curiously, the Gutenberg-Richter law for earthquakes potencies/moments, 
interpreted as statistics of dissipated energy, gives a very close value of the expo- 
nent e ~ 1.55 — 5/3 [28, 173]. This is not too surprising, however, in view of the 
essentially two dimensional nature of the fault friction/plasticity. 
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Figure 25: Log-log plots of the probability density for (a) dissipated energy, (b) strain increments, 
(c) avalanche durations. 



Due to very small scale of stress fluctuations the dissipated energy during an 
avalanche is approximately proportional to the size of the total plastic slip before 
the system stabilization. The probability distribution P{A) of the plastic strain 
increments 

A = N-'J2 f \i,Jdt, (43) 

is shown in Fig. 25(b) and, as expected, we again find the power law distribution 
with exponent yU = 1.6 ± 0.05. In this graph the abscissa is normalized by the 
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system size, and therefore one can see that the biggest events traverse the whole 
system. On the other hand, the lower bound in our scaling region corresponds to 
the smallest plastic strain increment which scales with the unit cell size (viscous 
scale is even smaller). This suggests that the cut off have purely geometric nature 
which is a signature of a scale free regime. 

In addition to energy it is also of interest to study the statistics of the avalanche 
durations [153, 109, 154]. The results shown in Fig. 25(c) suggest a power law 
scaling with the exponent r = 2.0 ± 0.1, see Fig. 25(c). This value can be 
compared with measurements in LiF micro-crystals where the duration exponents 
were found to be consistently bigger than 2 [154]. 

We now turn to the study of the power spectrum for dissipated energy E{t) 
which carries information about temporal correlations not only between avalanches 
but also inside a single avalanche [155]. The log-log plot of the computed power 
spectrum is shown in Fig. 26. It has a characteristic structure of the 1/ p noise 
with 7] = 1. Such signals, also known as pink or flicker noise, are encountered in 
many self organized critical systems ranging from astronomy to physiology [156], 
moreover the self organized criticality was originally proposed as an 'explanation 
of 1/f noise' [174]. 

It is important to point out that the time series obtained in our dynamical model 
result from slightly overlapping avalanches. This finding is consistent with the 
results of the numerical simulations in "running" sand pile models reported in 
[171], where the appearance of 1// noise was explicitly linked to the avalanches 
overlap at finite driving forces and where a very close RG estimate of the exponent 
?7 = 5/6 in 2D case was obtained. Our results also agree with the scaling relation 
r] = 3 — T derived in [167] under the assumption of random superposition of 
individual avalanches. 

Another important signature of a universality class is the average shape of 
avalanches with a given duration T. The avalanche shapes have been studied in 
different physical systems exhibiting critical behavior [157, 158, 49], in particu- 
lar, plastic avalanches are known to be slightly asymmetric (skewed) [159]. The 
average shapes of avalanches in our dynamic model are shown in Fig. 27. We 
observe that the avalanche shapes in the 2D model show much more variability 
with respect to durations than in the ID model. For instance, at longer durations 
the avalanche in our model are getting more symmetric, which is also a prediction 
of the mean field model [49], however, the agreement is not complete because 
some asymmetry always remains. In our calculations, we observed that the num- 
ber of units close to the marginal stability limit is increasing as the avalanche 
evolves towards its maximum which may be the origin of the tail responsible for 
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Figure 26: Power spectrum of the dissipated energy in the dynamic model is an example of a 1// 
noise. 
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Figure 27: Average avalanche shapes in dynamic model normaUzed by the magnitude and the 
duration. 

the observed asymmetry. In the dynamic model we were not able to perform the 
scaling collapse because of the large computational cost associated with explicit 
resolution of the fast events. This will be done later on in the automaton model. 

5.5. Spatial correlations in the dynamic model 

The spatial counterpart of the observed time correlations is the fractal struc- 
ture of dislocational patterns. Fig. 28 shows the energy density map f{9,$,) of 
the system after the fifth cycle. We observe that the energy distribution exhibits 
spatially separated local peaks associated with dislocation-rich regions which are 
rather stable in time. The distribution of these peaks is only partially guided by 
the quenched disorder and is largely a result of self organization. 

In order to characterize quantitatively the spatial complexity of this distribu- 
tion we constructed a level set representation for the energy density and identified 
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Figure 28: Spatial energy density distribution f{d, ^) after the fifth cycle. 



the boundaries of the spatially separated local peaks by selecting an irrelevant 
threshold (see the insert in Fig. 29). We then computed the energy density associ- 
ated with these regions Eq = ^ J2n /(^) obtained its probability distribu- 
tion which is a power law -P(-En) ~ E^'^ with exponent C = 1.45 , see Fig. 29. 
The obtained scaling is very robust spreading over five decades. 

Another way to quantify the fractal clustering is to compute the correlation 
function of the actual dislocation distribution shown in Fig. 30. We define 

C(r) = 2 (44) 
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Figure 29: The log-log plot of the probability density for the energies of dislocation rich regions. 
The insert shows the level set representation of the energy landscape shown in Fig. 28 




Figure 30: A fragment of the dislocation distribution in the dynamic model after the sixth cycle. 
Red and blue dots correspond to dislocations with positive and negative Burgers vectors, respec- 
tively. 

where N is the total number of dislocations and o/f^ is the number of pairs of 
dislocations with separation less than r [160]. If C(r) ~ with non-integer 
exponent D, this exponent is called the fractal dimension. In particular, for a 
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randomly distributed (not correlated) set of points, the fractal dimension D is 
equal to the dimension of space (D = 2 m our case). 

In our numerical simulations we observed that during the first loading cycle 
D ~ 2.0, which is expected given the random nature of the quenched disorder. 
With cycling, however, the long range correlations developed progressively and 
in the shakedown regime we recorded D ~ 1.74 independently of the initial dis- 
order, see Fig. 31(a). Interestingly, the dislocation patterns with D ^ 1.64 — 1.79 




(a) (b) 

Figure 31: (a) Log-log plot of the correlation function C(r); (b) Log-log plot of the box-counting 
measure N(s). 

have been observed experimentally in crystals with multiple slip systems; in sim- 
ulations with a single slip system fractal patterning has been previously linked to 
the possibility of dislocation multiplication [161] which is operative in our model. 

Another well-known method of identifying the fractal structure of a set is com- 
puting its box-counting dimension D [162] . It requires a covering of a set by 
boxes of size s and counting the asymptotics of the number of boxes N{s) as 
s — )■ 0. If N{s) = s^^ then D is the desired fractional dimension. Applica- 
tion of this method to the contour plot shown in the insert in Fig. 29 gives again 
D f» 1.75, see Fig. 31(b) which confirms the consistency of different measures of 
spatial correlations. 

In summary, our 2D dynamical model shows a broad interval of self organized 
critical behavior characterized by a series of critical exponents revealing temporal 
and spatial correlations (rj = 1, a = 1.6, D = 1.75 and ^ = 1.45). In the domain 
of parameters where we studied our system, these exponents do not depend on K 
and and are not affected by the level of disorder a. The system self-tunes to crit- 
icality by developing a highly correlated structure of persisting inhomogeneities 
which ensures effective energy dissipation. 
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5.6. The automaton model 

Despite the conceptual transparency of the dynamical model, the mechanism 
of reaching the critical regime remains obscure in this largely computational for- 
mulation. In an attempt to obtain a mathematically more transparent version of the 
model we reduce in this subsection our dynamic equations to an automaton model 
and compare the results of conceptually similar numerical experiments. As in ID 
case, we utilize the fact that in the quasi-static limit z/ — )■ viscous relaxation 
is instantaneous and the driven system can be viewed as almost always residing 
in the state of mechanical equilibrium. This dramatically simplifies the tensorial 
problem because instead of dynamic equations of visco-elasticity one needs to 
deal only with equilibrium equations of elasto- statics d^/du = 0. 

In order to solve the equilibrium equations analytically we again replace the 
smooth periodic potential shown in Fig. 12 by its piece-wise quadratic approxi- 
mation defined in each period [{d - 1 /2) (2^°) , {d+1/2) (2^°)] by 

Here d = 0, ±1, ±2, is the integer-valued spin variable describing a quantized 
slip. In the dynamic version of the model we we always implicitly assuming that 

e = 0.5. 

Suppose first that the lattice field dij is given. Then the equilibrium equations 
can be written in the form 

K[ui+ij + Ui^ij - 2uij] + [uij+i + Uij_i - 2uij - 2^°{dij - rfij-i)] 

This linear problem for the displacement field can be mapped into the Fourier 
space giving an algebraic equation 

[Ks+(q)C(q) + 5+(q)s;(q)]^(q) - 2r5;(q)rf(q) - ^(q) = 0. (46) 

where q = (g^., Qy) = {2ni/N, 211] /N) and all other notations have been already 
introduced in 5.2. By solving (46) we obtain 

2e%-(q)ci(q) + g(q) 

«(q) = — - — , (47) 

where A(q) = 2K(cos(gx) — 1) + '5~(q)s+(q). The Fourier image of shear strain 
can be then calculated as ^(q) = s~(q)u(q). 
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It is now straightforward to reformulate the dynamical problem as an inte- 
ger valued automaton. We first define the residual shear strain due to quenched 
disorder as ^/i(q) = Sy{q)(H{q)/\{q)). Then we observe that the variable 

AeM = a,-(^ + [a(q)]J)> 

representing shear strain associated with the slip must be confined between the 
limits 

-(d,, - l/2){2e) - [ihMi]+t < AC.M < {d., + l/2)(2e°) - [U^)]7]+t. 

One can see that the stability of a particular slip configuration dij is controlled by 
two random thresholds which also evolve with the loading (with time). 

When one of the thresholds is reached in one lattice point, the integer field dij 
is updated 

dij dij + Mij{d), 

where 

f 1, if Ae,,(t/) > -id,, - l/2)(2e°) - [a(q)] J + t, 
M,,(rf) =<;-!, if Ae,,(d) < [d,, + l/2)(2e°) - [Uq)]rj + 1 
otherwise. 

After each increment of loading At one has to check the stability of all units and 
continue the update until all units are stabilized. 

The update of the "slope" field A^j can be represented in Fourier space as 

Ae(q)^ Ae(q) + L(q)M(q). 

Here the kernel 

^(q) = . 2 '"'y.^, for q ^ 0, 

Ksm (gx/2) + sm {qy/2) 

can be viewed as the analog of the toppling matrix in the sand-pile models. The 
inverse Fourier image of this kernel in the real space L{x, y) is highly anisotropic, 
long-range and conservative (see Fig. 32). The conservativeness is understood 
here in the sand pile terms and means that after each re-distribution of shear strain 
(discharge), the total deformation ■ A^ij, which is controlled by the boundary 
conditions, remains unchanged. 
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Figure 32: (a) The real space representation and (b) the discrete map of the kernel (toppling matrix) 

L{x,y). 



In Fig. 32(b) we show the discrete map of the real space of the Green's func- 
tion L{x,y). The discrete map indicates the distribution of the strain for each 
discrete unit (spring) in response to an elementary slip which can also be viewed 
as a localized force dipole. From this map one can see see that if an element 
reaches a threshold an strain increase takes place in this unit. Away from this 
element, in the x direction, the strain is decreased in two neighboring elements 
but then it again increases reaching zero at infinity. In the y direction, the elastic 
strain simply decays. Observe that while the structure of our is different from 
the quadrupolar structure of the 2D elastic kernel in isotropic elasticity [15], the 
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main anti-ferromagnetic structure in the direction of the shp is preserved. In 
the limit K ^ oo the kernel I/(q) vanishes everywhere except q^. = 0, where 
1/(0, Qy) = 2^^ for Qy 7^ 0. One can see that the 2D model reduces in this limit to 
the ID model with the kernel (23) (NN version with /3 = 0). 

5.7. Numerical results for the automaton model 

As we have already noticed, the automaton representation greatly reduces the 
complexity of numerical computations which, in particular, allows one to deal 
with much bigger domains. In the automaton model fast depinning events are re- 
placed by jumps and outside the jumps the system evolves through a succession 
of equilibrium states. One can expect that such simplification of the model pre- 
serves the basic correlations between the avalanches but may affect the avalanche 
shapes. It may also affect the power spectrum particularly at high frequencies. In 
this section we begin a systematic comparison of the two models. 

The stress-strain curves in the automaton model are shown in Fig. 33. We 




Figure 33: (a) Strain-stress curve for the automaton model in the first cycle; the insert shows stress 
fluctuations, (b) Strain-stress hysteresis during the first 5 cycles. 

observe that the macroscopic response remains basically the same as in the dy- 
namical model. The small difference in the configuration of the hysteresis loops 
is due to the fact that in both models we compute the stress-strain relation only at 
some selected time steps. For instance, in the automaton model we know exactly 
the silent intervals and therefore can capture the strain/stress curve with higher res- 
olution. Instead, in the dynamic model we randomly chose our test points which 
leads to stronger fluctuations. 
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Figure 34: Distribution of dislocation cores obtained in the automaton model; red and blue colors 
correspond to dislocations of different signs. 

A typical steady dislocation pattern in the automaton model is shown in Fig. 34. 
One can see that the dislocations again form clusters without an obvious scale. A 
computation of the correlation function for this set of points gives a slightly bigger 
fractal dimension D ^ 1.9 ± 0.1 than in the dynamic model. 

The study of the temporal correlations in the automaton model can be based 
on the time series representing the dissipated energy (12). The typical signal in 
the shakedown state is shown in Fig. 35. Since now we have access to data for 
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Figure 35: Energy dissipation in the automaton model. 
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different domain sizes we can perform the finite-size scaling (FSS) analysis of 
the avalanche distribution [163]. In the critical state the power law probability 
distribution is expected to have a size dependent cut-off 

PiE) = E-^cp (^^^ , (48) 

where is a universal shape function and Ec is the cut-off energy with the scaling 
Ec ~ A^*^. According to the FSS hypothesis both, the exponents (6 and e), and 
the shape function Lp, define the universality class. If the FSS hypothesis is valid 




(a) (b) 



Figure 36: (a) Log-log plot of (T{q) versus q. The slope in the linear regime is 1.2. (b) Scaling 
collapse of the dissipated energy in the automaton model for different system sizes N. 

we must obtain for large E that < E'' >= f E'iP{E)dE ~ N''^'i\ where a{q) = 
5{q — e + 1). The exponent a{q) can be computed from the slope of the log- 
log plot of < -E^ > vs. and the slope of a{q) with respect to q is the cut-off 
exponent S. As we see from Fig. 39, the linear behavior starts from q ^ 0.6 
with a slope 5 = 1.2 ± 0.1. The exponent 5 ~ 1.2 ± 0.1 is close to the value 
5 ~ 1 obtained for plastic strain increments in [53, 169]. Once 5 is calculated, 
using < E >~ N'^^'^\ one finds the scaling relation (2 — e)6 = a{l) that leads 
to e ^ 1.6. The data collapse in the coordinates E/N^ and P(E)N^'^, presented 
in Fig. 36(b), shows an excellent agreement with the criticality hypothesis. The 
knowledge of the exponent 6 allows one to introduce a correlation length and study 
the size effect associated with collective interaction of dislocations [170, 56]. 

As we have already seen, one can also access temporal correlations through 
the analysis of the power spectrum for the dissipated energy. The power spectrum 
for our automaton is shown in Fig. 5.7. Observe that the 2D power spectrum is 



51 




Figure 37: Power spectrum in the automaton model has a power law structure P{f) / with 
r] = 1.57. 

different from the ID spectrum shown in Fig. 4.7. It is, however, similar to power 
spectrum one obtained in our dynamic model although with a larger exponent 
7] ~ 1.57 (which is, interestingly, very close to the value 1.59 ± 0.05 obtained 
for the 2D BTW sand pile [159]). The scaling part of the power spectrum is 
limited on the low frequency side by the duration of the longest avalanche. The 
absence of correlations in the low frequency region shows that the power spectrum 
reflects correlations within a single avalanches while missing correlations between 
different avalanches (see also Ref. [147]). 

We now move to the study of the avalanche shapes. Recall that in a critical 
state one can expect to see avalanches of all sizes and the average avalanche shapes 
must scale with their durations and magnitudes in a universal way [164]. This 
scaling hypothesis can be thoroughly tested in the automaton model. 

The first step is to average the avalanches of a given duration T and plot the 
resulting function < E(t)T > against the rescaled time t/T. In the critical state 
one should see the following scaling [165] 

< E{t)T >= T"'-^e{t/T), (49) 

where e is a universal scaling function and 7 is the exponent in the power law 
relation linking the average avalanche size < E{T) >=< //^^ E(t)dt > with its 
duration T 

< E{T) >~ TT. (50) 

The relation (50) is tested in Fig. 38, which confirms the power law scaling and 
gives the value of exponent 7 ^ 1.3. In Fig. 39(a), we present the average 
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Figure 38: Log-log plot of < E{T) > vs T. 



avalanche shapes E{t) for different durations before data collapse. The data col- 
lapse based on the relation (50) is shown in Fig. 39(b). The slightly asymmetric 
close to parabolic shape of avalanches is in agreement with the results of DDD 
simulations [159]. Notice that we presented only avalanches with durations in the 
limited range because below T = 20 the average avalanche shapes deviate con- 
siderably from the 'parabolic' pattern, while above T = 50 the data become noisy 
and the scaling collapse is less satisfactory. 




Figure 39: Average avalanche shapes of a given duration in the automaton model: (a) before and 
(b) after the scaling collapse. 

Another way to reveal the scaling of the of avalanche shapes is to compare 
avalanches with a given total dissipation Eq. In this case, one can anticipate the 
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following scaling relation 



< E{t)E, >= El~'^~e{t/E^). (51) 

where e is again a universal cut off function. The relation (50) suggests that the 
total dissipated energy scales with the avalanche duration as T ~ E^, where 
a = 1/7 ^ 0.76. 

In Figs. 40(a) and 40(b) we present our data for the average avalanche shapes 
at a given total energy before and after the scaling collapse. One can see that the 
scaling hypothesis works rather well. Interestingly, our average shapes (at a given 
total dissipation E^) are qualitatively similar to the average earthquake profiles at 
a given total moment [157]. The skewed shape is a manifestation of the fact that 
avalanches with a given total dissipation can have different durations. One can 
except that the maximum amplitude of the shorter avalanches is higher than of 
the longer avalanches with the same total energy. Then in the average shape the 
shorter avalanches will mostly contribute to the peak at small time while the long 
avalanches will be mostly responsible for the tail at large times which leads to the 
observed asymmetry. We observe the same asymmetry in the mean field theory 




(a) (b) 



Figure 40: (a) Average avalanche shapes at a given dissipated energy i?o before (a) and after (b) 
the scaling collapse. 

of depinning where the universal scaling functions e and e are known exactly 
for some classes of periodic potentials, for instance, in the piece-wise quadratic 
case e(x) ~ x(l — x) and e(x) ~ xe~^ [88, 166], with the energy duration 
scaling T ~ E^^"^ [88]. A direct comparison shows that neither the exponent a 
nor the average avalanche shapes match the predictions of the mean field theory. 
Similar analysis of the scaling functions for avalanches with given durations and 



54 



magnitudes in the cases of earthquakes and magnets, which also questions the 
validity of the mean field theory, can be found in [175]. 

Finally, one can check the consistency of the obtained data by varifuing the 
universal relations between the exponents e, r and 7. Thus, given that the proba- 
bility distributions P{E) and P(T) are related through P{T)dT ~ P{E)dE we 
can use the scaling relation < E{T) >~ T'^ to obtain 7(1 — r) = 1 — e[167]. 
This scaling relation works fairly well for our model showing the consistency of 
the obtained numerical data. 

6. Conclusions 

In this paper we developed a highly idealized model of crystal plasticity allow- 
ing one to study cooperative effects in dislocation dynamics. Quite remarkably, 
our toy model, based on a simple system of coupled ODEs, has shown an excellent 
agreement with experimental observations of fluctuations during plastic yielding 
and with previous much more elaborate numerical studies of avalanche statistics. 
The model confirmed robust temporal intermittency and spatial fractality of the 
shakedown regime and provided the first evidence towards the possibility of a 
flicker noise which remains to be checked in experiment. 

To make the occurrence of scaling in our model more transparent we further 
simplified the formulation by reducing the continuous dynamics to an integer au- 
tomaton. More precisely, we replaced the fast dissipative stages of dynamics by 
jump discontinuities and then linearized the elastic behavior inside individual en- 
ergy wells. The legitimacy of such model reduction has been proved rigorously 
only in ID case where the fluctuations have a Gaussian structure [101], and the 
generalization to higher dimensions appear to be nontrivial in view of the observed 
long range correlations. The derived automaton model, however, is much lighter 
numerically than the dynamic model which allowed us to perform the scaling col- 
lapse and to reveal the size effect (correlation length divergence) associated with 
collective interaction of dislocations. 

Based on the computed value of the critical exponents one has a temptation to 
conjecture that both models, the one with continuous dynamics and smooth poten- 
tial, and the one with discrete dynamics and piece- wise quadratic potential belong 
to the same universality class. This conclusion would be markedly different from 
the prediction of the mean field theory where smooth and cusped potentials lead to 
different universality classes [108, 168]. In fact, the dynamic and the automaton 
models show different exponents in the power spectrum. Also the fractal dimen- 
sions of the dislocation pattern in the two models are slightly different. Given that 



55 



the automaton model presents a drastic simplification of the discrete model it may 
be expected to capture only a limited interval of temporal and spatial scales. Dif- 
ferent behavior of the two models outside this interval may affect the maximum 
likelihood fit of the exponent and this may be the source of the observed disagree- 
ment. For instance, the remaining disagreement may be due to the small but finite 
overlap of the avalanches in the dynamic model while they are perfectly separated 
in the automaton model. More detailed studies are needed to resolve this issue 
definitively. 

An important remaining question is whether the toppling rules in the automa- 
ton are Abelian meaning that the outcome of the instability in multiple sites does 
not depend on the toppling order. We compared numerically the conventional 
updating strategies such as simultaneous updates of all unstable units, updates 
involving only the unit with the smallest value of strain or updates of randomly 
chosen units. We found that the microscopic configuration shows some small 
dependence on the choice of the strategy while the macroscopic manifestations, 
including the shakedown hysteresis loop and the statistics of avalanches (critical 
exponents) remain unaffected. This means that the model exhibits a weak (sta- 
tistical) form of Abelian symmetry. Another important property of the model is 
that the energy is lowered after each avalanche, which is explicit in the dynamic 
formulation and can be inferred for the automaton formulation. These two fea- 
tures, additional internal symmetry and the dissipative structure, are crucial for 
the possibility of the mathematical study of the origin of scale free behavior. 

To make our model more realistic and useful in mechanical engineering ap- 
plications it is of interest to develop a tensorial formulation allowing one to sim- 
ulate the critical behavior in the configurations with arbitrary geometry and with 
generic loadings. In particular, the complex inhomogeneous setting will give rise 
to interesting manifestations of the size effect associated with the growth of the 
correlation length from the microscopic values to the scale of the whole struc- 
ture. On materials science side, it would be of interest to include the effects of 
micro-cracking during cyclic loading associated with extreme dislocational pile- 
ups. Such coupling of plasticity and damage will open the way to the study of the 
crossover from plastic to 'spinodal' criticality associated with the formation of the 
macroscopic crack and the ultimate failure of a solid. 

More generally, the obtained reduction of a complex multi-particle dynam- 
ics to an integer automation suggests that similar development can be attempted 
for other distributed systems exhibiting criticality. The revealed inherent tempo- 
ral discretness can be a tool bringing dramatic simplification and transparency to 
otherwise impenetrable problems ranging from turbulent flows to the rheology of 
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cytoskeleton. 
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